%% This file allows us to solve a PDE (partial differential equation in 2D
% This example shows as we should proceed to solve a Partial Differential
% Equation (PDE) by using Local Max-Entropy (LME) shape functions.
% 
% The problem statement is:
%              | -Lap(f) = -2*y         in Dom,
%              | f = 0                  on the frontier of Dom,
% where Lap means Laplacian, and Dom is the domain, Dom=[0,1]x[0,1]
cd ..
addpath(path,'appHeatEq');

clear all
clear all
%clc
close all

%% PDE : Heat Equation (CODINA'S problem)
% Nomenclature
%   L2_@: L2 norm error
%   nn_@: nearest neighbors node list
%   n@  : number of samples or nodes
%   x_@ : samples or nodes coordinates
%   h   : nodal spacing
%   p_@ : LME shape functions
%   dp_@: LME shape function derivatives


L       = [1 1];
Nx      = 11;
Ny      = 11;
centre  = [0.5*L(1); 0.5*L(2)];
x_nodes = UniformGrid2D(Nx, Ny, L(1), L(2), centre);
nPts    = size(x_nodes,1);
h       = 1/(Nx-1);


%% Boundary conditions
% bc_homog = 0 homogeneus (the PDE is modified and appears a source term), 
% bc_homog = 1 non homogeneus (some nodes have BC imposed)
bc_homog = 1;
     

%% Gauss-Legendre Numerical Integration
% Numerical Integration
optGauss.orderGL = 4;     %cubature order: 2 (3 GPts), 4 (6 GPts), 6 (12 GPts)
% Sample points where the shape functions are computed
[x_gauss w_gauss conectivity gPts]= MakeGLSamples2D(x_nodes, optGauss);

figure(1)
plot(x_nodes(:,1),x_nodes(:,2),'ro',x_gauss(:,1),x_gauss(:,2),'bx')
xlabel('X')
ylabel('Y')
legend('node points','Gauss points')
axis equal


%% Computation of the shape functions and gradients
% Definition of locality parameters
gamma  = 1.6;
beta   = gamma/h/h;
beta_n = ones(1,nPts)*beta;

% Options for the computation of the shape functions
optLME.dim     = 2;           % spatial dimension
optLME.verb    = 1;           % information and plots: 0=OFF, 1=ON
optLME.grad    = 1;           % Computation of the Gradient 0:OFF 1:ON
optLME.hess    = 0;           % Computation of the Hessian  0:OFF 1:ON
optLME.TolNR   = 1.e-12;      % Newton-Raphson Tolerance
optLME.Tol0    = 1.e-6;       % Target Zero Tolerance
optLME.beta_n  = beta_n;      % locality parameters
range_n        = zeros(nPts,1);
for i=1:nPts
  range_n(i)  = max(2*h,sqrt(-log(optLME.Tol0)./beta_n(i)));
end
optLME.range_n = range_n;     % range definition

% adjacency structure with the nearest neighbors nodes to each sample point
% nodal shape parameter
%s_near = SamplesAdjacency(x_nodes,x_gauss,range_n);

s_near = samplesAdjacency_fem (x_gauss,conectivity,gPts);



% Local-max entropy basis functions computation
 optLME.s_near = s_near;
 %outLME = wrapper_lme(x_nodes,x_gauss,optLME);
 %p_samp  = outLME.p_samp;
 %dp_samp = outLME.dp_samp;


 % FEM basis function computation 
 outFEM = wrapper_fem (  x_nodes , x_gauss , conectivity, gPts);
 p_samp  = outFEM.p_samp;
 dp_samp = outFEM.dp_samp;
 

%% The numerical solution is computed
%id_bd: indices for the nodes on the boundary
nPts    = length(x_nodes);
ids     = (1:nPts);
id_bd   = [ids(abs(x_nodes(:,1))<1.e-6) ids(abs(x_nodes(:,2))<1.e-6) ...
           ids(abs(x_nodes(:,1)-L(1))<1.e-6) ids(abs(x_nodes(:,2)-L(2))<1.e-6)];
id_bd   = unique(id_bd);
u_nodes = ex23_SolveSystem(bc_homog,nPts,id_bd,x_nodes,s_near,p_samp,dp_samp,x_gauss,w_gauss);


%% Error in L2 norm
% compute value of u in the Gauss points
u_gauss = SamplesSolution(u_nodes,s_near,p_samp);

% analytical value of u in the Gauss points
u_sol   = ex23_AnalyticalSolution(bc_homog,x_gauss);

% L_2 norm of the error
L2_GL   = sqrt(sum((u_gauss-u_sol).^2.*w_gauss));

fprintf('\tNorm L2 of the error : %10.4e\n', L2_GL);

figure(2)
plot3(x_nodes(:,1),x_nodes(:,2),u_nodes,'ro',...
    x_gauss(:,1),x_gauss(:,2),u_gauss, 'bx',...
    x_gauss(:,1),x_gauss(:,2),u_sol, 'k+')
xlabel('X')
ylabel('Y')
zlabel('u(x)')
legend('node points','Gauss points num','Gauss points real')
